function [x,P] = rouwenhorst(mu, rho, sig, N)

% discretizes an AR(1) process using Rouwenhorst method
% - see Kopecky and Suen (RED 2010) for details
%
% Model:    x(t+1) = mu + rho x(t) + sqrt(sig) e(t+1)
% 
% mu is the intercept in AR(1) - zero if demeaned
% rho is the AR(1) coefficient/autocorrelation
% sig is Variance of residual

    if abs(rho)>=1; error('Process is explosive'); end

p = (1 + rho)/2;

psi = sqrt((N-1)*sig/(1-rho^2));

y = linspace(-psi,psi, N);

theta{2} = [p 1-p; 1-p p];

if N>2

    for n=3:N

        theta{n} = p*[theta{n-1} zeros(n-1,1); zeros(1, n)] ...
            + (1-p)*[ zeros(n-1,1) theta{n-1}; zeros(1, n)] ...
            + (1-p)*[ zeros(1, n); [ theta{n-1} zeros(n-1, 1)]] ...
            +    p *[ zeros(1, n); [ zeros(n-1, 1) theta{n-1}]];
        theta{n} = theta{n}.*[ones(1,n); .5*ones(n-2,n); ones(1,n)];
    end

end

P = theta{N};

x = y' + mu/(1-rho);